Introduction

The Urban Heat Island (UHI) refers to the phenomenon where urban areas experience higher temperatures compared to their surrounding rural areas. This effect is primarily caused by the built environment of cities, which absorbs and retains more heat, as well as the greater presence of human heat sources such as cars and air conditioners (Rizwan et al., 2008). In the city of Utrecht, the UHI is evident through a temperature gradient, where the highest temperatures are observed in the city centre, gradually decreasing towards the rural areas (Brandsma & Wolters, 2012).

Elevated temperatures resulting from the UHI can have an impact on plant communities. Increased temperatures can increase the abiotic stress plants experience, resulting in an alteration between competition and facilitation dynamics within plant communities (Olsen et al., 2016). Consequently, the interactions between dominant and subordinate plant species may shift from facilitative to competitive. Subordinates may experience reduced survival due to this increased competition. These dynamics can ultimately lead to a reduction in plant diversity (Olsen et al., 2016). Given that the UHI phenomenon raises local temperatures, a negative effect of the UHI on plant diversity can be expected.

In our research, we investigated this impact in the city of Utrecht by examining road verges in different parts of the city. We assessed the plant species richness in flower-rich hay meadows of different UHI. Our hypothesis is that plant diversity declines with increasing intensity of the UHI. We expect sites with higher UHI to have less plant diversity than sites with lower UHI.

Methods

We collected data in plots of flower-rich hay meadows located in different areas of Utrecht, ranging from the city centre to the eastern regions and near Utrecht Science Park. All of the sites were roadside verges that are managed by the municipality. Our selection of the specific sites was based on the mowing policy implemented by the Gemeente Utrecht (Maaikaart, 2023). We specifically used meadows classified as Bloemrijk hooiland 1x maaien (GR3+G12), as they represent fairly natural and accessible grasslands. By choosing sites with similar grassland types and mowing policy, we aim to eliminate any potential confounding variables associated with variations in grassland type or mowing frequency. To increase our statistical reliability, we selected 12 different sites for replication purposes. Figure 1 shows the location of the plots through Utrecht.


Figure 1: Plot locations in the city of Utrecht. The locations are numbered to ensure that single data points can be traced back to individual plots. Figure from https://goo.gl/maps/umGwFhXZe3ewnn3T8
Figure 1: Plot locations in the city of Utrecht. The locations are numbered to ensure that single data points can be traced back to individual plots. Figure from https://goo.gl/maps/umGwFhXZe3ewnn3T8


The number of plant species were determined using morphological differences. Specifically, we visually identified and counted the number of species present within a 1x1 meter plot at each site. To reduce biases, the plots within each site were carefully selected to be representative of the overall plant diversity at the site. For each plot, we determined the intensity of the UHI at that location, defined as the difference in temperature on that site compared to rural area’s. For this, we used the report of Brandsma and Wolters (2012) on the UHI effect in Utrecht. This report contains a figure of the modelled UHI effect in the city and its surroundings (Figure 1). Using the colour scale, we assessed the UHI for each site, incorporating the mean of each UHI class in our data.


Figure 2: Spatial distribution of the mean night-time UHI intensity relative to the background temperature of the surrounding rural areas. The black curve represents the transect along which the measurements were taken. Retrieved from Brandsma & Wolters (2012).
Figure 2: Spatial distribution of the mean night-time UHI intensity relative to the background temperature of the surrounding rural areas. The black curve represents the transect along which the measurements were taken. Retrieved from Brandsma & Wolters (2012).


Results - statistical analysis


The data set


Data import

The data is given in the file “Data UHI”. We import it into a data frame which we call “Data”.

Data <- read.csv("Data UHI.csv",  sep=";")
Data
Data structure

We examine the structure of the data set.

str(Data)
## 'data.frame':    12 obs. of  3 variables:
##  $ Plot       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ UHI        : num  1.45 1.15 0.8 0.35 1.05 1.75 1.65 1.85 1.75 0.8 ...
##  $ Diversiteit: int  19 10 17 13 12 6 9 6 12 11 ...

The data set contains three vectors: “Plot” specifies the plot number; “UHI”, specifies the intensity of the urban heat island and “Diversiteit” contains count data of the total number of plant species of each plot.

We calculate some descriptive statistics.

summary(Data)[,2:3]
##       UHI          Diversiteit   
##  Min.   :0.3500   Min.   : 6.00  
##  1st Qu.:0.9125   1st Qu.: 9.75  
##  Median :1.3000   Median :12.00  
##  Mean   :1.2500   Mean   :12.00  
##  3rd Qu.:1.6750   3rd Qu.:14.25  
##  Max.   :1.8500   Max.   :19.00
sd(Data$Diversiteit) ## Standard deviation of Diversiteit
## [1] 3.977208

We can also summarize these statistics in a boxplot:

library(ggplot2)  ## Load package which contains the function ggplot

ggplot(Data, aes(factor(0), Diversiteit)) +
  geom_boxplot(width = 0.4) +
  geom_jitter(width = 0.06) +  ## Add jittered data points
  labs(title = "Plant diversity",
       y = expression("Amount of plant species" ~ ("m"^-2)),
       x = "") + 
   theme_bw() +
  scale_x_discrete(breaks = seq(0, 2, by = 0.2),
                     labels = function(x) ifelse(x == 0, "", x)) +
  scale_y_continuous(breaks = seq(6, 19, by = 2),
                     labels = seq(6, 19, by = 2))

Data visualization

We now plot Diversiteit along UHI.

ggplot(Data, aes(x = UHI, y = Diversiteit)) +
  geom_point() +
  labs(title = "Plant diversity along UHI",
       y = expression("Amount of plant species" ~ ("m"^-2)),
       x = "Urban Heat Island (UHI)") +
  theme_bw() +
  scale_x_continuous(breaks = seq(0, 2, by = 0.2),
                     labels = seq(0, 2, by = 0.2)) +
  scale_y_continuous(breaks = seq(6, 19, by = 4),
                     labels = seq(6, 19, by = 4))


Analysis

Assumptions

We perform a general linear model framework and test whether the assumptions of linear regression are met. This will be done by the means of diagnostic plots.

Result_lm <- lm(Diversiteit ~ UHI, data=Data)

par(mfrow = c(1,2)) # Make a 1 x 2 panel graph
plot(Result_lm, 
 which = 1:2) # Only plot the first two graphs

The first plot checks for the assumption of linearity and homoscedasticity. The red line is not reasonably straight and horizontal, therefore the assumption of linearity is violated. The assumption of homoscedasticity, however, does seem to be respected, as the variance of the errors is fairly constant. The second plot checks for normality of the error distribution. The residuals in the second figure are roughly on the same line, therefore this assumption is respected.

We could use a generalized model to make the assumptions of homoscedasticity and normality of the residuals fit even better. We first try to transform the data to make the relationship more linear.

# Square root transformation
Result_lm_transf <- lm(Diversiteit ~ UHI + I(sqrt(UHI)), data = Data)

par(mfrow = c(1,2)) 
plot(Result_lm_transf, 
 which = 1:2)

# Log transformation
Result_lm_transf <- lm(Diversiteit ~ UHI + I(log(UHI)), data = Data)

par(mfrow = c(1,2)) 
plot(Result_lm_transf, 
 which = 1:2)

# Perform a linear regression on the transformed data
Result_lm_transf <- lm(Diversiteit ~ UHI + I(UHI^2), data = Data)

par(mfrow = c(1,2)) 
plot(Result_lm_transf, 
 which = 1:2)

The transformations do not improve the linearity of the relationship. However, because we have very few data points, we think we can proceed with our analysis without transforming the data. We can check this with the Box-Tidwell test, which tests the linearity assumption and examines the significance of the interaction term. We use the car package to perform the Box-Tidwell test:

library(car)
## Loading required package: carData
boxTidwell(Diversiteit ~ UHI, data=Data)
##  MLE of lambda Score Statistic (t) Pr(>|t|)
##         11.314             -1.2679   0.2366
## 
## iterations =  18

The non-significant result of the Box-Tidwell test suggests that there is no significant evidence to reject the null hypothesis of linearity. In other words, there is no strong indication of a violation of the linearity assumption in the relationship between the predictor variable (UHI) and the response variable (Diversiteit). We therefore consider the assumption to be satisfied.

We will now use a generalized linear model with Poisson error distribution and a log link function, which is best suited for count data. We account for overdispersion by using a quasi-poisson regression.

Result_glm <- glm(Diversiteit ~ UHI, family= quasipoisson(link = 'log'), data=Data)

summary <- summary(Result_glm)
summary
## 
## Call:
## glm(formula = Diversiteit ~ UHI, family = quasipoisson(link = "log"), 
##     data = Data)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.8650     0.2570  11.148 5.82e-07 ***
## UHI          -0.3121     0.2021  -1.544    0.154    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.253032)
## 
##     Null deviance: 15.022  on 11  degrees of freedom
## Residual deviance: 12.068  on 10  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

We extract the slope of the regression with 95% confidence interval and the p-value.

summary$coefficients[2,1] ## Slope of the regression
## [1] -0.3121464
confint(Result_glm)[2,] ## Confidence interval for the regression slope
##       2.5 %      97.5 % 
## -0.70648491  0.08709861
summary$coefficients[2,4] ## Pr(>|t) for the regression slope
## [1] 0.1535674

We plot Diversiteit along UHI, adding a regression slope with 95% confidence interval that corresponds to our glm model.

ggplot(Data, aes(x = UHI, y = Diversiteit)) +
  geom_point() +
  geom_smooth(method=glm,  method.args=list(family = quasipoisson(link = 'log'))) +
  labs(title = "Plant diversity along UHI",
       y = expression("Amount of plant species" ~ ("m"^-2)),
       x = "Urban Heat Island (UHI)") +
  theme_bw() +
  scale_x_continuous(breaks = seq(0, 2, by = 0.2),
                     labels = seq(0, 2, by = 0.2)) +
  scale_y_continuous(breaks = seq(6, 19, by = 4),
                     labels = seq(6, 19, by = 4))


Conclusion

We can conclude that there is no significant relationship at the 95% confidence level between the number of plant species per plot and intensity of the Urban Heat Island (slope and 95% CIs = -0.31 (-0.71 - 0.087, p = 0.154). However, our results do show a trend towards lower species richness with increasing UHI. This is what we would have expected if increasing temperatures do change the interactions between plant species from being facilitative to competitive. Future research can collect more data to find a significant result. We still cannot say, however, that there is a causal relationship between plant species richness and UHI, as we conducted an observational study and UHI was not directly manipulated. Other factors, like nitrogen emissions, water availability and other soil characteristics would likely also differ between different parts of the city of Utrecht and have an effect on the measured differences in plant diversity.


References

Brandsma, T., & Wolters, D. (2012). Measurement and Statistical Modeling of the Urban Heat Island of the City of Utrecht (the Netherlands). Journal of Applied Meteorology and Climatology, 51(6), 1046–1060. https://doi.org/10.1175/JAMC-D-11-0206.1

Maaikaart. (03-2023). Retrieved 1 June 2023, from https://gu-geo.maps.arcgis.com/apps/webappviewer/index.html?id=7373671535e8426f8945113a9ee8e97d

Olsen, S. L., Töpper, J. P., Skarpaas, O., Vandvik, V., & Klanderud, K. (2016). From facilitation to competition: Temperature-driven shift in dominant plant interactions affects population dynamics in seminatural grasslands. Global Change Biology, 22(5), 1915–1926. https://doi.org/10.1111/gcb.13241

Rizwan, A. M., Dennis, L. Y. C., & Liu, C. (2008). A review on the generation, determination and mitigation of Urban Heat Island. Journal of Environmental Sciences, 20(1), 120–128. https://doi.org/10.1016/S1001-0742(08)60019-4